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SUMMARY 

Three-dimensional unsteady transonic flow through an axial turbomachine stage is 
described in terms of a pair of two-dimensional formulations pertaining to orthogonal 
surfaces, namely, a blade -to -blade surface and a hub-to-casing surface. The resulting 
systems of nonlinear, inviscid, compressible equations of motion are solved by an expli- 
cit finite -difference technique. Separate computer programs have been constructed for 
each formulation. The blade -to -blade program includes the periodic interaction between 
rotor and stator blade rows. Treatment of the boundary conditions and of the blade slip- 
stream motion by a characteristic type procedure is discussed in detail. Harmonic anal- 
ysis of the acoustic far field produced by the blade row interaction, including an arbitrary 
initial transient, is outlined in an appendix. Results from the blade -to -blade program are 
compared with experimental measurements of the rotating pressure field at the tip of a 
high-speed fan. The hub-to-casing program determines circumferentially averaged flow 
properties on a meridional plane. Blade row interactions are neglected in this formula- 
tion, but the force distributions over the entire blade surface for both the rotor and stator 
are obtained. Results from the hub -to -casing program are compared with a relaxation 
method solution for a subsonic rotor. Results are also presented for a quiet fan stage 
designed by the National Aeronautics and ^ace Administration, which includes transonic 
flow in both the rotor and stator and a normal shock in the stator. 

INTRODUCTION 

Flow through a high-speed fan or compressor is highly three-dimensional and can 
include complex shock wave systems. In addition, flow through a complete stage consist- 
ing of a rotor and stator or a fan preceded by inlet guide vanes is unsteady even in the 
rotatii^; frame of reference. Effects of viscosity and turbulence are known to be signifi- 
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cant; in fact, the turbulent wakes may be the predominant source of aerod 3 mamic inter- 
action between the rotating and stationary blade-rows and responsible for the associated 
noise. Vortex filaments are known to stream from the rotor tips and undoubtedly inter- 
act with the wall boundary layer. Sufficiently far downstream, turbulence generated by 
the first blade row may encompass the entire flow. . C^culation of transonic flow field 
solutions in high-speed turbomachinery stages is clearly one of the most formidable 
challei^es to present-day capabilities in computational aerodynamics. 

A traditional approach to solution of. this complex problem has been, taken. 
Description of the inviscid flow field is addressed first, boundary layer and turbulence 
effects being. svg)erimposed as small perturbations. The basic system of equations which 
is solved numerically consists of the complete nonlinear equations of motion for an invis- 
cid compressible gas: 

+ V . pV = 0 (1) 


where 
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Dt P . 
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t time ■ 

V velocity vector 
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p density 

The only essential simplifications introduced pertain to reduction of the spatial 
dimensions of the problem from 3 to 2. ’ Solution of the full three-dimensional problem is 
reduced to a pair of two-dimensional solutions on orthogonal surfaces. A separate pro- 
gram tailored to the particular aspects of each formulation has been developed. Both 
computer programs which will be discussed herein are intended for analysis of a com- 
plete stage, that is, rotor and stator, although an isolated blade row can also be treated. 

A harmonic type analysis of the acoustic far field due to blade row interactions has 
also been developed. This analysis provides a direct coupling between the numerical 
near -field solution and the acoustic far field, with a particular view toward character iz-' 
ii^ the acoustic far field. This aspect of the flow model will be only briefly outlined 
herein, but a complete description is given in reference 1. The boundary -layer and wake 
representations utilize standard integral methods and it is assumed that the boundary 
layer and wake are quasi-steady. The reader is again referred to reference 1 for a com- 
plete description. 

SYMBOLS 

An,Bn,Cn,. . . Fourier coefficients 
a speed of soimd 

b stream sheet thickness or blade thickness 

t centerline 

c chord 

E total internal energy 

e internal energy 

H total enthalpy 

h enthalpy 

M Mach number 
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m 


meridional distance along a blade -to -blade stream surface 
Nj number of blades in ith row 

n blade -to -blade distance 




n 


P,Q,R boundary conditions at interface of numerical near -field and acoustic far- 
field solutions 


p pressure 

R gas constant 

r radial distance from axis of rotation 

S entropy 

s distance aloi^ the stream path i 

T temperature 

t time , , 

U a reference velocity ,, , y 

u meridional component of velocity; ^ 

u velocity component parallel to slipstream or blade surface , 

V velocity vector 

V circumferential component of velocity 

V velocity component normal to slipstream or blade surface. 

Vr radial component of velocity; also circumferential yelocity comp.onent in 

rotating coordinate system 

vz 
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axial velocity component 



Vg circumferential velocity component 

X meridional distance along a blade-to-blade stream surface 

X . distance along slipstream or blade surface 

y . circumferential distance 

y .distance normal to slipstream or b^de surface 

z axial coordinate 

^n>^n acoustic propagation coefficient 


y ratio of specific heats 

0 circumferential angle , 

di,92 blade surface coordinates 

p density 

(p angle between slipstream or blade surface and meridional plane 

ij/ ratio of mass flow to total, mass flow 

angular speed of rotor 
w vorticity; also frequency 

Subscripts: 

o initial or reference condition 
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00 


free stream 



A circumflex (^) over a symbol denotes in average value. A bar over a symbol 
denotes vector quantities except in appendix where it denotes time to frequency transform. 

i . . , : . 

BLADE -TO-BLADE ANALYSIS 


In the blade -to -blade formulation the previously stated equations of motion are 
e^ressed in a curvilinear coordinate system alined on an axisymmetric stream surface 
as shown in figure 1. This stream surface Is considered to have small but finite thick- 
ness b(z) and variable mean radius from the axis of rotation r(z). The velocity com- 
ponent normal to the stream surface is neglected; as a result, a two-dimensional approx- 
imation to the flow field is produced. (See refs. 2 and 3.) When the. m,0 coordinates 
are transformed to a rotating system by 

X = m (7) 

y = r(0 - fit) (x = Constant) (8) 

and the circumferential velocity component is transformed by ... 

Vj. = V - fir (9) 


then the following system of governing equations result: 
^ + ±(pu) + -^(pvr) ^ 

at dx^ ' ay^ rb dx 


( 10 ) 


^ * -^(pu2 . p) + ^uvr) . . 


ay' 

ay' 


_ _puf drb 
rb dx 


— + — (puVr) + — (pVr^ + p) = - - pufi . 

at ax'^ ay'^ ^ rb dx ^ ‘ 


dx 


* #KHr) = ^ - pua . V,) 

where the relative total energy and total enthalpy are defined by 
Hj. = H - firv 


( 11 ) 

( 12 ) 

(13) 

(14) 


Et 


Hr -K 


(15) 


The terms oh the right-hand side of the equations result from the variations in cross- 
sectional area and radius of the stream surface, which are intended to account for the 
effects of variations in hub and casing radii. The terms on the left-hand side of the 
equations correspond to the familiar. set of two-dimensional unsteady equations of motion 
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of an inviscid compressible gas. Note that the relative total enthalpy is not the same as 
^thp total enthalpy defined on the basis of velocity components in the rotating frame. 
Numerical solution of this system of equations at interior grid points is accomplished by 
the MacCormack algorithm (ref. 4). Systematic rotation of the order in which the non- 
centered differences are evaluated is employed to minimize any bias in the solution due 
to the ^ternating directions of the noncentered differences. No artificial damping or 
stabilization is used. 

■ As shown in figure 1, the computational domain is divided into a maximum of seven 
se^entsj- not all of which need be included in every case. The grid network extends 
axially from an inlet station to a discharge station and circumferentially across one 
blade -to -blade passage. The inlet station is located either 1 axial chord length iqjstream 
of the first blade row, in which case domain 1 is deleted,, or an arbitrary distance 
upstream exceedii^ 1 chord length. Placement of the discharge station can be selected 
in the same manner as that employed for the inlet station. Domains 2 to 6 are each lin- 
early mapped into a unit square which is spanned by a rectangular grid network. In 
domains 1 and 7 a linear stretching of the axial coordinate is used to map these domains 
into unit squares. The axial grid spacing in domains 1, 4, and 7 is determined by the 
locations of the axial bovmdaries of . these domains. The lateral boundaries of domains 1 
and 2 lie on extensions of the mean camber line. The lateral boundaries of domains 3 . . 
and 5 lie on the blade surfaces. The instantaneous locations of the blade slipstreams 
form the bovmdaries of domains 4, 6, and 7. Domains 1 to 4 are attached to the first 
blade row and domains 5 to 7 are attached to the second blade row either of which may be 
selected as the rotating row; It is assumed that the number of blades in the second row 
equals or exceeds the number in the first row. 

It is emphasized that a periodic solution due to the aerodynamic interaction of the 
rotating and stationary blades is anticipated. The formulation is thus a numerical coun- 
terpart of the problem for which Kemp and Sears (ref. 5) obtained an analytic solution 
pertaining to thin, slightly cambered blades of low solidity in an incompressible flow. 
Linearized solutions have more recently been obtained for compressible flows, but the 
avithors are not aware of any other attempt to develop nonlinear solutions for, the periodic 
blade -row interaction problem at transonic or supersonic conditions. 

In connection with the periodicity of the subject problem, two main points of depar- 
ture from other numerical solutions for transonic airfoils or cascades should be recog- 
nized. First, the slipstreams are moving surfaces of discontinuity across which jumps 
in tangential component velocity and in total pressure can occur. Only static pressure 
and the component of velocity normal to the surface must be continuous. It should be 
noted that the jumps in tangential velocity across the slipstreams are related to the 
unsteady variations in lift of the blades and therefore cannot be obtained from the con- 
servation form of the equations of motion through a limiting process as‘ in the case of 
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shock waves.' -Treatment of the slipstreams as surfaces of discontinuity in the present ,i. 
model is therefore warranted for two reasons: It allows, attainment of an accurate, peri 7 ; 
odic solution without requiring a very large grid point density to approximate the slip- 
stream discontinuities, and tracking of the slipstreams is necessary to determine the 
trajectories of the viscous wakes which diffuse outward from the inviscid slipstreams. 

The second point is that in the case of an unequal number of blades in the two rows, the ’ 
angular period of the circumferential variations in the flow field is not the width of a 
blade -to -blade passage but it is the circumference divided by the difference between the 
number of blades in the rotor and stator. Furthermore, the flow pattern rotates :with an 
an^lar velocity which, in general, is a multiple of the wheel speed. Numerical repre-^ 
sentation of this periodicity condition pertaining to the lateral boundary points of the grid 
network, as well as to those points along the interface between domains 4 and 5, is ' ' 
accomplished by a cyclic procedure which is discussed later. . - 

BOUNDARY CONDITIONS 

The calculation of boundary points and indeed the exposition of proper boundary 
conditions is facilitated by recasting the equations in the form of characteristic compati- , 
bility relations pertaining to a quasi one -dimensional unsteady wave system (as suggested 
by Mbretti and Abbett (ref. 6 ) and Serra (ref, 7)), . Although actual numerical implementa- 
tion of the characteristic formulation is far more complex than the finite -difference pro- 
cedure used in the interior points and can possess certain drawbacks, such as Inconsis- 
tency with the interior point solution, it is nevertheless adopted here for the particular 
advantages it offers with respect to the slipstream, inlet, and discharge point calculations. 

Consider first the inlet station sketched in figure 2, As is. well known, the naonien- 
tum and energy equations can be rewritten in the form: . 

^ = 0 or S = Constant on ^ ^ = M (18) 

Dt ; u Vr • . ' 

= 0 or Constant on ^ = At (17) 

Dt . ^ u Vr 

By assuming that no flow reversal occurs at the inlet station, the values of the entropy 
and ratio of vorticity to density at the inlet station are, therefore, solely properties of the 
incoming flow and may be specified a priori; both the entropy and vorticity of the Incom- 
ing flow are assumed to be zero. Three options have been considered with regard to 
physical interpretation of the inlet boundary conditions. First, if the Inlet station repre- 
sents an open end of a finite length duct, the static pressure can be specified as 

P = P_oo 
^94 



Second^ if> it is an arbitrary station in a duct of infinite length, the Riemann' invariant on 
the downstream traveling wave can be specified as 


2a 

V - 1 


.+. u 


_ 2a -00 

r - 1 


+ Uoo 


(18b) 




Use of this, relationship implies that the outward traveling waves are one-dimens^nal 
(i.e., either-planes or helices). As a third option the numerical solution can be matched 
to the acoustic far -field analysis at the inlet station, in which case all acoustic modes 
are properly accoimted for. This procedure is outlined in an appendix. In any of these 
cases, the. solution of the inlet boundary points is completed by first using a compatibility 
re,lation bn the, pp stream runnii^ wave which originates within the computational domain 
at point C in figure 2; that is, 


I A log p - Au = 




9y 


a 9 log P 
y dy 


- a ^ - iia 

ay dx 


At 


^ori M = u - a and ^ = Vo) (19) 

and second by solving the following angular momentum equation, which only involves 
spatial derivatives of the dependent variables in the circumferential direction, by the 
MacCormack finite-difference algorithm: , . 


(20) 

This. characteristic compatibility relation may be interpreted as pertaining to the projec- 
tion of the true characteristic conoid on a reference plane which is alined normal to the 
inlet station and translates in the circumferential direction with an arbitrary velocity Vq, 
for example, Vf from the previous time step. The acoustic far -field analysis can be 
used to. replace equation (18a) or (18b). 


8Vr’ 

at - ay 
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+ Vr 


avf 
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uS2^ 

dx 


A similar procedure is employed at the discharge boundary points but here the 
entropy and vorticity are determined by tracing a particle path from within the computa- 
tional domain (point B in fig. 2) to the boundary point. The system of equations pertain- 
ing to a discharge boundary point are stated below for the case of either a finite length 
duct or ah iiifinite duct 

P^Poo ' , . ^21a) 


2a 2a«, 

- u = 22^ - Uoo 

y 1 y - 1 


(21b) 
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Use of the acoustic model would replace equation (21a) or . (21b) 
S = S* 


0) _ 


P* 




ion ^ = ^ = At) 


(22) 


S A log p + Au 


- Vo)(|j + ^ gy 
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a and 


au av». 1 ap „ w 

ay ^ ay P ay ^ P* 


* + ufi^ 
dx 


^ = v 
At 
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(23) 

(24) 


The boundary condition at the blade surface points is simply vanishing of the com- 
ponent of velocity normal to the surface. At the trailing edge the Kutta condition is sat- 
isfied by requiring the pressure to be continuous and the velocity component normal to 
the mean of the camber line and slipstream to be zero. On the slipstream the pressure 
and normal component of the velocity must.be continuous. .Implementation of these condi- 
tions is facilitated by recasting the equations into a characteristic form similar to that. ^ 
described; however, in this case the reference planes are normal to the surface and 
translate in the streamwise direction, as shown in figure 3, Combination of the continu- 
ity equation and normal momentum equation results in the following pair of compatibility 
relations: 


a A Iqg p ± y Avr ;= (aQj ± yQ 2 )At 


on 


|Z=Vr±.a. and |S = Uo: 


(25) 


where 


Ql = - 


(n - uo)U2£P + is + yu 
^ dx dx dx 


(26) 


Qo = - 


(u - Uq)-^ - uO cos <^> ^ 
dx dx 


(27) 


The energy equation is stated as 

DS = o 

Dt 


(on Ax = u At and Ay = v At) 


S = Constant 


(28) 
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and the streamwise momentum equation is solved in the LaGrangian form' 

Dq _ 9 log P o dr 

Dt " y 9s " q dx 


(29) 


where q^ =\u2 + and ds/q = dx/u = dy/v = dt. 


It. is pointed out that the form of the compatibility relations given by equation (25) 
provides an algebraic solution for the quantities which are continuous across the blade 
slipstream, namely, the pressure and normal velocity. Thus, the boundary conditions on 
the blade slipstream can be satisfied without iteration, other than that necessary to locate 
the characteristic geometrically (points A and D in fig. 3) by successive approximations, 
of which two are usually sufficient. However,, at the trailing ec^e. an iteration is required 

• • ii • ■ • 

to determine the slipstream angle which satisfies the Kutta condition. 


The overall scheme for imposing the boundary conditions along the blade surface 
and slipstream points is shown schematically in figure 4; The time axis projects verti-:. 
cally out of the page in this, figure. The dashed lines represent the intersections of the- 
translating reference planes with the axisymmetric stream .surface and .the intersections ‘ 
of the:particle paths with the stream surface during a-, time step At, < • 


PERIODICITY CONDITION 


Illustration of the nature of the cyclic procedure devised to enforce the periodicity 
of the solution can best be accomplished with respect to the following simplified configu-’ 
ration. Consider first a stage having three rotor blades and three stator blades. This ' . 
configuration is shown in figure 5 in both axial and cascade projections. At time /to. all” 
rotor and stator blades are alined, whereas at time to + At the rotor has moved through 
a fraction of a revolution, and none of the blades are now, alined. It is cle^ in this case 
that the geometric conditions which determine the flow through the stage are identical in 
each blade -to -blade passage at any time.?- In this case .the solution along an exterior grid 
row 6 can be equated to that along the interior grid line )3 and similarly that along 
exterior line a can be equated to that along interior line y at any instant. Consider 
now the case with three blades in the stator and four blades in the rotor as shown in fig- 
ure 6. At time to rotor blade 2 is alined with stator blade b, whereas at time to + At 
rotor blade 3 is alined with blade c. In this case the geometric conditions pertaining to 
the passage between blades a and b are obviously different from those for the passage 
between blades b and c at any time. However, it may be noted thatjthose pertaining 
to passage be at' to + At are precisely the same as those which pertain to passage ab 

?^It is assumed, of course, that the boundary conditions imposed at the inlet and dis- 
charge stations are spatially uniform so that the blade geometry provides the only scale 
for circumferential variations. 
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at the previous time to. Therefore, the flow, conditions along exterior grid line 6 at 
time to + At can be equated to those along interior grid line /3 at the earlier time t©. 
However, in this case those along exterior grid line a at time to + At cannobbe 
equated to those occurring in passage ab at time to, but must be equated to_those 
occurring along line y at an earlier time. Thus a. phase shift is introduced in applica- " 
tlon of the. lateral boundary conditions. The necessary boundary information /is acquired , 
during the passage of time, and therefore the desired periodicity is attained asymptoti- 
cally in time. , 

, NUMERICAL EXAMPLE: BLADE -TO-BLADE PROGRAM 


Results from the blade -to -blade program have been compared with data for a high- ^ 
speed. (1 500 -fps) fan tip section for which experimental data are reported in reference 8. 
The casing wall was instrumented with an array of fast response pressure gages, from, 
which a contour plot of the rotating pressure field around the tip section was recon- 
structed. This fan was preceded by a set of guide vanes and followed by a row of stators. 
However, the unsteady interaction was neglected in this case and only the rotor was con-- ^ 
sidered. The grid network was very coarse and consisted of 9 grid rows in the circum- ' 
ferential direction and li in the axial direction in each of 3 domains, that is, a total of 
297 grid points. 


The e3q>erimental pressure contour plot is reproduced on the left-hand side of fig- 
ure 7. The data indicate the presence of an oblique shock off the leading edge of the 
upper blade which reflects off the lower blade and reimpinges on the upper blade near the 
trailing edge. A lambda'(\) type shock is apparently formed on the aft portion of the 
\ipper blade because of the boundary-layer separation. The isobars constructed from the 
numerical solutions are shown on the right-hand side of this figure. The numerical ' " ^ 
results exhibit qualitatively similar behavior, although boundary -layer effects have not 
been included and the grid is admittedly very coarse. Although a qualitative correlation 
is apparent with respect to the main features of the flow field, a quantitative comparison 
is somewhat difficult. Therefore, the experimental isobars have been used to construct 
pressure distributions along the suction and compression surfaces of the blade and along 
a mid -channel line. The accuracy of the data obtained in this manner may be somewhat 
suspect, but the agreement between the data and the numerical solution shown in figure 8 
is considered to be very encouraging. 

Numerical solutions for interacting blade rows have thus far been limited to ideal- 
ized test case configura.tions for which no comparisons with other solutions or experi- 
mental data are available. However, resjxlts for the full stage consisting of fan and sta- 
tor tip sections from reference 8 are included in reference 1. 



HUB-Tb-CASING ANALYSIS 


•i- 


Attention is now shifted to the second program, which considers a hub-to-casing ' . 

, • • ... » ^ 
stream surface. The coordinate system and grid network is illustrated in figure 9. The 

finite -difference grid lies on a meridional plane extending from hub to casing and from an' 
inlet station to a discharge station. Circumferential variations are removed in this case ' 
by integration of the governing equations with respect to d from one blade to the next 
and defining average properties over this ang^ar interval. The angular velocity compo- 
nent and angular momentum equation are retained because of the presence of a pressure 
force exerted by the blades. The effects of blockage of the flow area due to the blade 
thickness and boundary -layer displacement thickness are included. Multiple blade rows 
can be considered but the effects of periodic interactions between the blade rows are 
necessarily neglected because of integration of the equations with respect to the angular' '• 
variable. 

In this analysis the basic system of equations given by equations (1) to (6) are stated 
In cylindrical coordinates, multiplied by d0 and integrated from 0j(r,z) to d 2 {r,z), 
which are the surface coordinates of two adjacent blades, as indicated in figure 9. Out^ 
side a blade row, the integration interval is taken as 2n/Ni. The blade -to-blade passage 
width is defined by 
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(Within a blade row) 
(Outside a blade row) 


(30) 


Circumferentially averaged values of the dependent variables are defined by 

■*02 


= n'^Nj C ” pr d0 

J01 


(31) 


PVz 


= n -%1 

Jfl. 


• prv2 d0 


(32) 


It is assumed that the blades are thin and sharp and hence the Tocar blade surface 
angles can be replaced by the mean camber line angles: 


9^1 ££2 _ M 

3z 9z 9z Camber line 


(33) 


1 ~ 2 80 

9r 8 r ®*‘ICamber line 

Furthermore, differences between root -mean -square (rms) values and mean sqiiared 
values are neglected: 

- Vz^l « 


VrVz - VrVz « a‘ 


(35) 

(36) 


The following system of equations is thereby obtained^: 
a(pn) 


^ + -2_(pnvz) + —ipnvr) = 0 
at 8z ^ ar'^ 


(37) 


+ 4 (pnvz 2 ) + ^(p„vr) = -n ||> 4 p(r ||) 


(38) 


■l^pnvr) + A(pnvrv^) * -^priyr^) . -n g + Ap(r M) 


(39) 


^(pnrv 0 ) + ^(pnrvgVz) +^(pnTVgVr) = -r Ap (40) 

^(pnE) + ^(puVzH) + ^(puvrH) = -rO Ap (41) 

where , 

Ap(r,z) = p(r,z, 02 ) - p(r,z,0i) (42) 

Thus the variable n represents the circumferential distance around the annulus, reduced 
by the cumulative blockage due to all blades, at fixed r and z. The term Ap repre- 
sents the cumulative pressure differential across the blades, that is, the pressure differ- 
ential across each blade times the number of blades. The terms involving Ap in the 
momentum equations represent the three components of the pressure force exerted by the 
blades. The term on the right-hand side of the energy equation is the work performed by 
the rotating blade row. 


2 >rije superscript notation to denote average quantities is dropped in the following 
discussion. 
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Within a blade row the velocity vector is assumed to remain tangent to the mean 
camber surface at all times: 

“ V0 = Or + vz ^r II) + vr(r ||) (43) , 

The tangency condition is used within the blade rows to determine the angular velocity 
component, and the blade pressure differential Ap is obtained from the angular momen- 
tum equation. Outside the blade rows Ap = 0 and the angular momentum equation is 
used to determine the ai^lar component of velocity, 

/ 

Representation of the inlet and discharge boundary conditions and numerical solu- 
tion of these boundary points follows the general approach described in connection with 
the blade -to -blade program. However, since a steady, rather than periodic, solution is 
sought in this case, somewhat less care need be taken in modeling a physically correct 
boundary condition. In particular, the "infinite duct" condition discussed previously has 
been replaced in this case by prescription of the total pressure at the inlet station. Sim- 
ilarly, prescription of the components vorticity of the incoming flow can be replaced by 
direct statement of the flow aisles or of the radial and circumferential velocity compo- 
nents. At the discharge boundary the static pressure is specified. Solution at the bound- 
ary points along the hub and casing surfaces is accomplished by restating equations (37) 
to (41) in a body-oriented coordinate system similar to. that used to derive equations (25) 
to (29). In this case, the boundary condition V = 0 replaces the normal momentum 
equation. The same finite -difference procedure as used at the interior points is used to 
accomplish the solution of the remaining members of the system of equations. Noncen- 
tered differences are used for the derivatives normal to the walls, which only involve 
gradients of the normal velocity component. The accuracy of this procedure has been 
foxmd to compare very well with that of the interior point solution, with considerable less 
complexity than the characteristic procedure used in the blade -to -blade program. 

NUMERICAL EXAMPLES: HUB-TO-CASING PROGRAM 

The present hub-to-casing program is analogous in many respects to the relaxation 
method program developed by Katsanis and McNally at NASA Lewis Research Center. 
Their prograha MERIDL (ref. 9) solves the stream fimction equation on a meridional plane 
through a blade row for steady subsonic conditions by use of a finite -difference method. 
The comparison between the present program and their program has been carried out for 
a case in which their solution should be very accurate. The rotor configuration selected 
corresponds to the test case used by Katsanis and McNally in reference 9, 

The relative swirl angle is shown in figure 10. This angle is defined as that 
between the velocity vector and the projection of the velocity vector on a meridional plane, 
measured in a rotating frame of reference. The three curves pertain to the hub surface. 
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the tip or casing surface, and a mid -channel surface. Within the blade row the relative 
swirl angle along the hub and casing is completely determined by geometric constraints. 
However, everywhere else it is obtained from the solution for the three velocity compo- 
nents. The agreement between the two programs is' considered to be quite satisfactory. 

The. magnitude of the velocity vector in the rotating frame is shown in figure 11. 

The mid-channel values were deleted from this figure for clarity. The only significant 
difference between the two solutions occurs near the leading edge. Both programs allow 
grid columns to cross the leading and trailing edges in an arbitrary fashion. However, 
Katsanis and McNally's program accounts for the effects of bluntness pf the leading edge 
in some detail whereas the.present prc^ram assumes that these edges are sharp. 

. Next, a transonic fanjstage designed by NASA Lewis, Research Center ..for the. Quiet . ; 
Fan Program has been considered. Designated the QF-1 stage it combines an 1100-fps.;_ 
tip speed rotor with a stator having highly "leaned" blades (up to 45^ at the tip). A range 
of stator positions relative to the rotor location are possible with this stage; position VI 
was selected in this case. The rotor and stator are relatively close in this position, less 
than 1 chord length apart. : ' _ 

The pressure distributions along the hub and casing surface are shown in figure 12, ■ 
and the absolute Mach numbers are displayed in figure 13. The supersonic region' which 
develc^)S in the stator along the hub; due to the condbination of the hub curvature and blade 
thickness, is terminated by a normal shock. This shock is spread over about 4 or 5 grid 
points, or about 25 percent to 30 percent of the 16 grid points which cover the stator 
axially in this case. Only 10 grid points cover the rotor tip section; consequently, a 
shock in the rotor would be difficult to detect with the present grid point density, A rapid ' 
conipression is evident on the aft portion of the rotor tip section, which corresponds to a ’ 
reduction in relative Mach numbers from 1.12 near the leading edge to 0.64 at the trailing- 
edge. Therefore, the rotor -tip compression substantially exceeds the normal -shock 
compression by itself. , . , 

COMPUTER EXECUTION TIME AND. STORAGE REQUIREMENTS 

The blade -to -blade program will fit in small-core memory of a GDC 7600 computer, 
that is, about 160K octal words, with a maximum of 1000 grid points (not including the 
exterior points required for the boundary point calculations), exclusive of the storage 
required for the periodic boundary data needed for unequal numbers of blades. In its 
present form, disk storage is used for the periodic boundary data, although use of the 
large -core memory would undoubtedly be more efficient. The blade -to -blade program 
requires approximately 2 x 10-4 second per grid point per time step for execution on a 
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CDC 7600 computer by usir^ the FTN (opt = 2) compiler to generate thie binary code. The ‘ 
results presented in figures 7 and 8 required less than 1 minute of execution time. 

The hub-to-casing program also fits in„sniall-core memory with a maximuna of 
1600 grid points. It requires approximately 4 x 10"^ second per grid point per time 
step for execution. The additional time relative to the blade -to -blade program is 
believed to be associated with calculation of the blade pressure differential Ap, solution 
of an additional momentum. equation, and continuous reevaliiation of the maximum permis- 
sible time step. (A variable time step is not allowed in the blade -to -blade program 
because of the procedure for storing and retrieving boundary data.) The subsonic rotor 
case discusseil in connection with figures 10 and 11 required less than r rriimite of execu- 
tion time with a 27 by 17 grid network. The transonic stage results shown in figures 12 
and 13 required approximately 7 minutes of execution time. v 

CONCLUDING REMARKS 

A blade -to -blade formulation and a hub-to-casing formulation have been developed 
for analysis of transonic unsteady flow through an axial turbomachine stage and imple- 
mented in two computer prc^rams. Both employed an explicit finite -difference technique 
for solution at the interior grid points, and a characteristic type procedure at the inlet 
and discharge boundaries. The blade -to -blade prc^ram can treat periodic interactions 
between rotating and stationary blade rows, and particular attention has been given to 
correct representation of the blade slipstreams and their motion due to unsteady blade 
loading. A comparison with experimental measurements of the rotating pressure field of . 
a high-speed fan tip section is considered to be very encouraging. The computer execu- 
tion time for this case was very modest, less than 1 minute on a CDC 7600, and use of a 
higher grid point density to improve the numerical accuracy is therefore practical. 

The hub-to-casing program compares favorably with a relaxation time solution for 
a flow condition when the latter should be very accurate. ResiUts have also been obtained 
for a quiet fan stage which includes transonic flow in both the rotor and stator and a nor- 
mal shock in the stator. The program resolves the shock reasonably well, although a 
higher grid point density would probably be beneficial in this case. 

These programs offer a substantial improvement in the predictive c^abilities 
available to aerodynamicists involved in design and evaluation of high-speed turboma- 
chinery stages. However, it is clear that a complete description of the three-dimensional, 
imsteady,. turbulent flow prevailing in such stages will require continued development of 
the computational models. 
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APPENDIX 


ACOUSTIC FAR-FIELD ANALYSIS 

Under the conditions typically prevalent in highly loaded transonic fan or compres- 
sor stages, the linearized, small -perturbation approximations to the equations of motion 
cannot be expected to be descriptive of the flow in the vicinity of the blades. Thus, 
recourse is made to the numerical solution of the complete nonlinear system of equations 
as discussed in connection with the blade -to -blade program. However, sufficiently far 
from the blade rows, the amplitude of the flow disturbances will decay to acoustic levels 
and the linearized, small -perturbation approximations will be descriptive of the far field. 
Therefore, an intermediate region in which both analyses are valid should exist at some 
distance from the blades. The inlet and discharge stations of the blade -to -blade compu- 
tational domain can serve as the interfaces between the near-field (numerical) and fair - 
field (acoustic) analyses. The present far-field analysis is formulated with respect to an 
infinite duct model, namely, outgoing waves should propagate without reflection. It dif- 
fers, however, from conventional inlet duct analyses in that the signal may begin with an 
arbitrary, transient, associated with the deviation of the assumed -initial data in the near 
field from the periodic -solution which is sought as the asymptotic limit in time. There- , 
fore, the. acoustic analysis must recognize that a transient signal will occur during ^ 
startup and that a simple harmonic time dependence, which is the usual basis of inlet duct 
acoustics, cannot be assumed. The analysis should allow the transient to radiate, outward 
without reflection, and should be capable of identifying the attainment of a periodic solufri 
tion by the growth of discrete harmonic components in the solution. 

The inlet and discharge stations indicated in figure 1 as the axial boundaries of 
domains 1 and 7 (or possibly of domains 2 and 6 if 1 and 7 are deleted) form the axial ' , 
boimdaries of the acoustic far field. However, the lateral boundary should extend; over • 
that fraction of the circumference which corresponds to the fundamental period of the . 
stage configuration (that is, an integer number of blade r to -blade passages). This will -■ 
require storage and retrieval of numerical data along this interface in the same manner ; 
as is performed along the interface between domains 4 and 5 of the. near field. v , 

As discussed in connection with the characteristic procedure used at the inlet and’ 
discharge boundaries, the compatibility relation on the outward running, waves (eqs,. (19) 
and (23)) provide a connection between the combination of pressure and axial velocity on 
the boundary points at time t -h At and the known interior point solution at time t. This 
information provides the mechanism for matching the acoustic far-field solution with the 
numerical near-field solution on a point -by -point, step-by-step basis. Since the blades, 
are capable of producing a vorticity field which will convect downstream, additional infor- 
mation is necessary to define the downstream far field. (Recall that the inlet flow is 
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assumed to be irrotational.) K the standard small -disturbance approximations are 
employed, the numerical data to be provided at the inlet consist of 

P(y,t) = - (u - u_oo) (Al) 

P-00**’-oO 

and at the discharge station 

R(y,t) = v (A3) 

where the subscripts denote the reference states at x — ±», which are not necessarily 
the same. Equations (Al), (A2), and (A3) form, in effect, the boundary conditions for the 
far -field analysis, from which the instantaneous values of the pressure perturbation and 
velocity perturbation at the inlet and discharge stations are obtained. 

With these preliminaries in hand, attention is now focused on the acoustic far -field 
analysis. In the following discussion the subscript ( )q will be used to denote the refer- 
ence conditions for either boundary, that is, x — ±<», and all variables without subscript 
will refer to perturbations with respect to the reference state, namely, p = p - Pq and 
U = U - Uq. . 

Since the present effort is addressed toward a cascade formulation, the governing 
equations are written in a two-dimensional Cartesian coordinate system. It is noted that 
the two-dimensional problem could be described by a solution of the wave equation alone, 
were it not for the fact that the time dependent force distribution on the blades is capable 
of producing a convected vorticity field (as well as a corresponding entropy field, which, 
however, is not relevant to the present problem). It will be seen that the convected vor- 
ticity field (characterized by a solenoidal velocity component) does not contribute to the 
acoustic pressure field by itself, as distinguished from the irrotational component of the 
velocity field which is directly coupled with the acoustic pressure. Therefore, the veloc- 
ity perturbation field downstream of the discharge boundary station and upstream of the 
inlet boundary station is characterized by the sum of an irrotational velocity yector uj 
and a solenoidal velocity vector U 2 , which satisfy the linearized conservation equations^: 

^ + Po ^ o ^^'^=0 

at 

Pq ^ + Vp = 0 (A5) 

^Note that U 2 = 0 upstream of the inlet; however, the analysis is developed with 
respect to the more gener^ case pertaining to the flow downstream of the discharge 
station. 
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^ a2dP-PdS_n 

(jt -ao — - — — -0 


dt Cy dt 


(A6) 


where 


X U2 = ^ 0 

A = 3 + U-^ 
dt 9t ax 

• U2 = 0 

ui = Fui + fvj 

X UJ = 0 

U2 = Fu2 + fv2 

= ui + U2 

U = Moa.Q 


It can be implied from these equations that solutions for p and uj are purely radia- 
tive (propagating acoustically in a coordinate system convecting at the mean flow veloc- 
ity U) whereas solutions for U2 are purely convective. These results can be 
expressed in terms of a Fourier integral representation: 


p(x,y,tf 


PoaoAn(«^r 

/ui(x,y,t)' 


Bn(‘^) 1 

vi(x,y,t) 

n 

Cn(w) 


^i(cot-j3nx) ^ 
2ir 


(A7) 


fu2(x,y,tn V- 
lv2(x,y,t)J ^ 


2n 


(A8) 


The boundary data can be expressed as 

P(y,t) 

|Q(y,t) 


Pn°(<^)' 


R(y,t) 


^ = ^e’^“ny y** |QnV)S ei<^t d^ 


(A9) 




where a bar oyer: a.S3mibol denotes the time to frequency transform, and the svqperscript 
o indicates transformation from spatial location y to spatial harmonic n which will 
be discussed later. 


Substitution of these integral forms into the governing equations gives 
^n 


Cn = B„ 


(AlO) 
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(All) 


En = 


_ Ct> 




" oT- " 


(A12) 


An + ®n + ~ Qn° (A13) 

« 

(A14) 


The value of A^, Snd thus Bn, Cn, Dn, and En, can thereby be expressed in terms of 
the transform of the boundary data as ' 

An = - 5 ^ ‘ ' (^Qn-^UQtnRnO) ' (A15) 

co^ + (ao - + Uaooio^ ■ - . 

The propagation coefficient /3n is deterniined by substituting the integral relation for, 
pressure into the convected wave equation as 


^n = 


-cuMq ± 


u}^ - Q!n2ao2 (l - Mq^) 
ao(l - Mq2), 


1/2 


(A16) 


where the + sign refers to downstream propagation (the discharge boundary) and the 
- sign to upstream propagation (the inlet boundary). 

The selected representation of the solutions and boundary conditions as Fourier 
series in the y-direction is appropriate for enforcement of the periodicity boundary condi- 
tion pertaining to spatial variations in this direction. The coefficient an is defined 
accordingly as 



(A17) 


.where Y is the fundamental period of the stage cascade configuration. In addition, the 
fact that the boundary data are specified at a discrete number of grid points, say N, on 
the boundaries implies that the Fourier series can only include N terms, that is, 
n = 0, 1, 2, . . ., N - 1. Since the distance y to each point can be written as mY/N, 
where m = 0, 1, 2, . . ., N - 1 also, the Fourier series can be expressed in the st^d- 
ard Discrete Fourier Transform (DFT) flotation: ‘ ^ ' 


rp(y,t)" 

Q(y,t) > 

.R(y,t). 



^Pn°(t)' 

Qn°(t) 

J^n°(t) 


g-27Tinm/N 


(A18) 
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The inverse DFT is then 


Pn°(t^ 

N-1 

Pmft)' 

jQn°(t) 



Rn°(t) 

m=0 

Rm(t) 

s- J 


(A19) 


The desired solution for the pressure perturbation on the boundaries is aecom- 
plished after some manipulation as 


0(t) + J„0(t) * 

°(t) 

where P^W refers to the value at one of the N boundary grid points on the inlet or 
discharge boundary. Use of DFT techniques for Pm(t), Qm(t), and Rm(t) and similar 
DFT expansions for Hm, Jm, and then leads to the convolution: 



N-1 


V o- 

p^(t) = p^ao ) e-2^inm/N ) 

n=0 lKn°(t) * Pn 


r^n(t) * Qm-n(t) + •Tn(t) ♦ Pm-n(tr] 

p (t) = ^ y < 

^ lK„(t) *Pm-n(t) ;; 


N Lj 
n»0 


(A21) 


on the discharge and inlet boundaries, respectively. Thus a double convolution over both 
time and distance (in the y-direction) is required. The functions Pm(t)> Qm(^)> 

Rm(t), therefore, represent the point sources of time -varying strength which are alined 
along the considered boundaries, spatial resolution being consistent with the number of 
grid points specified. The functions Hui(t), Jni(t), and Km(t) are the duct response 
functions defined by 


"HmCtf 

N-1 i 


/ Jni(t) 


V(t) [ 


I- n=0 


Km(t) 


Kn°(t) 

J 


g -27Tinm/N = 


\ «(« 

1 



N-1 


1 ° 


V(t) i 


n=l 


i6(t) 

V. 


Kn°(t) 


-27rinm/N 



(A22) 


”n“<^) = 5T r (« - T JSJelfi’- fS 

.V. n2(l - Mq^) 27t 


Jn°(T) = - r (n - 


(A23) 


(A24) 
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Kn°(^) = -7 h r + MoTn)ei«^^ (A25) 

n(l -Mo2)(1 + Mo)‘^-«^ 

where 6 is the Dirac delta function and 

- (l 


i 

A groip of subroutines including an efficient Fast Fourier Transform routine has 
been developed for use in conjunction with the blade -to -blade computer program to carry 
out the indicated convolutions numerically. Results have been thus far limited to test 
cases with a simple harmonic input signal. For example, in one of the calculations the 
discharge station was assumed to be divided into 8 intervals covering a total circumfer- 
ential distance of 0.1 foot. The selected reference (average) Mach number of the dis- 
charge flow was 0.8. The flow was assumed to be ir rotational so that only one input 
function Q(y,t) was required, and the second input fxmctioh R(y,tj could be considered 
as a response function (that is, it was calculated from Q(y,t)). The input function 
Q(y,t) = cos^Ot - with Q(y,t) = 0 for t < 0, was selected for this case, where 

T = 2irNaot/Y, n = 1, 2, . . ., N, N =8, Y = 0,1 foot, ao = 103 fps, and n = 1. The 
input function Q and response fxmction R at n = 1 are plotted in figure (Al). It 
should be noted that the response function is initially out of phase with the input function 
because of the assumption that Q = 0 for t < 0. . However, the effect of the transient at 
T = 0 dies quickly, and after about 1/3 millisecond the response function closely approx- 
imates the input function and indicates the desired harmonic solution is being approached 
asymptotically. The nondimensional perturbation pressure is plotted in figure (A2). The 
complete history is shown for the point n = 1, whereas the history of the points n s 2 
and 3 is only shown at early times, where a difference in amplitude as well as phase 
exists. At later times (that is, after about 1/3 millisecond) the pressure solutions at 
the various grid points only differ noticeably by the phase angle corresponding to the 
input faction, as the differences in amplitude asymptotically decay. 


Tn = 




a 


T= « 


aot 


o = 2l 
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Figure 1. - Blade -to -blade coordinate system and grid network. 














Figure 6. - Illustration of cyclic algorithm for stajge with imequal number , 
of blades in stator and rotor. Nj = 3, N2 = 4. 


614 




Figure 7.- Comparison of experimental and numerical pressureydistributions 

for 1500 fps rotor. .. 
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Figure 13.- Variation of Mach number aloi^ hub and casing thrmigh NASA QF-1 stage. 
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